We simulate a random sample of 1000 Normal\((\mu=0, \sigma=5)\) IID random variables.

n <- 1000
x <- rnorm(n, mean = 0, sd = 5)

The log-likelihood function can be computed with the following function:

l <- function(x, mu, sigma){
  sum(dnorm(x, mean = mu, sd = sigma, log = TRUE))
}

We plot the log-likelihood function for \((\mu, \sigma)\) given \(\vec{x}\), saving the values in l_surface (code hidden).

library(plotly)
plot_ly(z = ~l_surface) %>% 
  add_surface(x = sigma, y = mu)

We compute the MLE's \((\widehat{\mu}_{\text{MLE}}, \widehat{\sigma}_{\text{MLE}})\) for \((\mu, \sigma)\):

xbar <- sum(x)/n
sigma_hat <- sqrt((1 / n) * sum((x - xbar) ^ 2))
c(xbar, sigma_hat)

## [1] 0.02755272 4.82191388

The MLE \((\widehat{\mu}_{\text{MLE}}, \widehat{\sigma}_{\text{MLE}})\) are the parameter values that best explain/fit the observed data x.